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Abstract. 

A two-dimensional fracture model where the interaction among elements is modeled 
by an anisotropic stress-transfer function is presented. The influence of anisotropy on 
the macroscopic properties of the samples is clarified, by interpolating between several 
limiting cases of load sharing. Furthermore, the critical stress and the distribution 
of failure avalanches are obtained numerically for different values of the anisotropy 
parameter a and as a function of the interaction exponent 7. From numerical results, 
one can certainly conclude that the anisotropy does not change the crossover point 
7 C = 2 in 2D. Hence, in the limit of infinite system size, the crossover value j c = 2 
between local and global load sharing is the same as the one obtained in the isotropic 
case. In the case of finite systems, however, for 7 < 2, the global load sharing behavior 
is approached very slowly. 
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1. Introduction 

For many years, the scientific community has shown great interest in the fracture of 
composites materials under imposed external stresses [fl |2[ [3] . By now several aspects 
of this process are well understood but a definite and complete physical description 
has not been made yet. Furthermore, the huge technological impact of composites 
materials has led to a continuous development of new models and theoretical approaches 
[H El E] . In fracture mechanics, the average mechanical properties of the specimen are 
commonly considered to be the input data for material modeling. [TJ, [2], [3] . Nevertheless, 
heterogeneous materials, such as fiber-reinforced composites, present widely distributed 
local mechanical properties. Thus, analytical approaches are very limited and numerical 
simulations have become an indispensable tool in this field. On the other hand, the latest 
developments in statistical mechanics have led to a deeper understanding of breakdown 
phenomena in heterogeneous systems [H [2], [3] . 

Fiber reinforced composite materials exhibit a large variability of ultimate 
macrosopic properties. Heterogeneity and anisotropy, in the micro- meso- and macro- 
structure of the composite, result in a complex scenario of damage mechanisms. 
Basically, the damage mechanisms include fiber breakage, matrix cracking and yielding, 
fiber-matrix debonding and delamination jH El El C] . In this framework, the simulation 
of the composite behavior may be achieved by the statistical modeling of the micro- 
structure and the development of the relation between micro-structure and macro- 
behavior d El El U\. 

Until now the modeling of the fracture of laminar composites has been based 
on finite element calculations FEM and some micro- mechanical models [H El [TO] . 
This modeling, has clarified that under pure shear loading the overall response of the 
sample is controlled mainly by the resin response. In fact, the scientific community 
recognizes that FEM techniques provide an excellent tool for predicting composite 
performance in controlled loading conditions. However, the continuous nature of FEM 
models usually makes them unable to describe the local damage evolution; which is the 
primary micromechanical process. Therefore, a local approach is mandatory for fully 
understanding this complicated process, from physics viewpoint. 

During the last two decades Monte Carlo simulations have been used to numerically 
study stress redistribution in 2D and 3D for different fiber arrangements [H El El CI- 
As a results, several aspects of composite fracture, when the external load is parallel to 
the direction of the fibers, have been clarified. Nevertheless, the fracture process and 
damage evolution of anisotropic systems, such as laminar composite materials subjected 
to shear external stress, are far from being well understood. 

A very useful approach to the fracture problem are the well-known Fiber Bundle 
Models (FBM), introduced long time ago by Daniels [UJ and Coleman [TJ] and subjects 
of intense research during the last years [lJ[TIl[lja[Tia[lIl[^ 

FBMs are constructed so that a set of fibers is arranged in parallel, with each one having 
a statistically distributed strength. The specimen is loaded parallel to the fiber direction 
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and the fibers break if the load acting on them exceeds their threshold value. Once the 
fibers begin to fail, several load transfer rules can be chosen. The complex evolution of 
internal damage and its associated stress redistribution are the most important factors 
to take into account in the accurate prediction of material strengths. The simplest case is 
to assume global load sharing GLS, which means that after each fiber failure, its load is 
equally redistributed among all the intact fibers remaining in the set. Otherwise, in local 
load sharing LLS the overload is only transferred to the nearest neighbors. This case 
represents short range interactions among the fibers. However, in actual heterogeneous 
materials stress redistribution should fall somewhere between LLS and GLS. 

In this paper, a generalized discrete model, where the interaction among elements 
is described by an anisotropic stress-transfer function, is introduced. By varying the 
anisotropy strength and the effective range of interaction we interpolate between several 
limiting cases of load sharing. The work is organized as follows. In section [2] the model 
and the way in which simulations are carried out are explained in detail. Numerical 
results are presented and discussed in section El The conclusions are given in the final 
section. 

2. The Model 

The fracture of fiber-reinforced composites is characterized by a highly localized 
concentration of stresses at initial cracks. Anisotropic laminar reinforcement prevents 
the nucleation of small cracks and the propagation of damage. In this way, the final 
collapse of small cracks in a critical cluster is avoided, retarding sample failure. 

In materials science, the Weibull distribution [2] has proved to be a good empirical 
statistical distribution for representing sample strengths, P{o~) = 1 — e~^~^ P . p is the 
so-called Weibull index, which controls the degree of threshold disorder in the system 
(the bigger the Weibull index, the narrower the range of threshold values), and a D is 
a reference load which acts as unity. On the other hand, in continuous homogeneous 
materials, the load profiles around a local damage area can be well fitted by a power 
law, 

^ ~ (1) 
where Act is the stress increase on a material element at a distance r from the crack 
tip. The above general relation covers the cases of global and local load sharing, 7 — » 0, 
and 7 — > 00, respectively. The transition from these limiting cases has been successfully 
described in isotropic systems [16j. In global load sharing approach, the strength of the 
sample can be computed analytically as vols = ^ = (pe) -1 / p for a Weibull distribution 
and o~gls — ^ — \ f° r a uniform distribution P(er) = 0~/a o , of breaking thresholds. 

Starting from these results [16] a discrete model with anisotropic load sharing 
is introduced. In the model, it is assumed a two-dimensional square-lattice of N 
elements, each one having a strength taken from a given cumulative distribution P(o~) 
and identified by an integer i, 1 < i < N. Thus, to each element i, a random threshold 
value o~i th is assigned. 
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The system is driven by quasi-statically increasing the load on each element. The 
element i breaks when its stress <7j is equal to its threshold value o~i th . Hence, the 
minimum value of o~k — o~k th in the set / of all unbroken elements, 

So-™ = min [a k - a k J (2) 

k e I 

defines the load increment Scr™ 1 ™. The quasi-statically load increasing is then performed 
by adding the amount of load 5cr™ n to all the intact elements in the system. Following 
this approach, all intact elements have a nonzero probability of being affected by the 
ongoing failure event, and the additional load received by an intact element i depends 
on its Axij and Ai/y from the element j which has just been broken. Furthermore, an 
anisotropic interaction is assumed between elements such that the load received by an 
element j, due to the failure of i, follows the relation: 

F(Axij, Ay ij} 7, a) = Z> (aA Xij 2 + (1 - a)A^/)" 7 , (3) 

where Ax^ and Ay^ are their relative distances on the x-axis and y-axis, respectively. 
7 and a are adjustable parameters and the value Zj is always given by the normalization 
condition, 

Zi=l/J2 {aAx id 2 + (1 - a)Ay l3 2 ) ^ , (4) 

The sum runs over the set I of all unbroken elements. We assume periodic boundary 
conditions, which means the largest value of Ax^ and Ayij is ^jp, where L is the linear 
size of the system. We note here that the assumption of periodic boundary conditions 
is for simplicity. In principle, an Ewald summation procedure would be more accurate. 
In Eq.(j3]) the extreme cases 7 — > and 7 — > 00 also correspond to global load sharing 
and local load sharing, respectively. Strictly speaking, for all a, the range of interaction 
covers the whole lattice. When changing the anisotropy factor a, one moves from a 
completely anisotropic case (either a = or a = 1) to the isotropic load redistribution 
(a = 0.5). 

Following this approach, a failing element transfers its load to the surviving elements 
of the set. This may provoke secondary fractures in the system which in turn induce 
tertiary ruptures and so on until the system fails or reaches an equilibrium state where 
the load on the intact elements is lower than their individual strength. In this latter 
case, the external force is increased again and the process is repeated until the material 
macroscopically fails. The size of an avalanche A is defined as the number of broken 
elements between two successive external drivings. Hence, during an avalanche of failure 
events, an intact element i receives the excess load from failing elements j at each time 
step. Consequently, its load increases by an amount, 

o-i(t + r)= ^{t + r-l)F{r ijn ,a), (5) 

jeB(r) 

where the sum runs over the set B(t) of elements that have failed at time step r. Thus, 

T 

o~i(t + T) = °"i(^o + T ) is the total load element i receives during an avalanche 

T=l 
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initiated at to an d finished at to + T. In this way, when an avalanche ends, the external 
load is increased again and another avalanche is initiated. The process is repeated until 
no intact elements remain in the system and the ultimate strength of the material a c , is 
defined as the maximum load the system can support before its complete breakdown. 

3. Simulation Results 

The mechanical properties of the bundle and the statistics of internal damage events were 
studied numerically by varying the anisotropy strength (a) and the range of interaction 
(7). Large scale numerical simulations in 2D were executed. Several system sizes 
(L = 65, 129, 257, 513, 1025) were considered and simulations were performed over at 
least five different realizations for the biggest system L = 1025 and five thousand for 
the smallest one L = 65. We recorded the avalanche size distribution D(A), and the 
ultimate strength of the samples which is always normalized by the characteristic 
value a of the cumulative threshold distribution P(er) = P(^-) ■ 

It is important to remark, that by using an isotropic power law stress redistribution 
Aa a dd ~ 1 a crossover point was observed [IB] . Hence, two distinct regions were 
distinguished, over the domain of 7. For small 7, is independent on L, which 
corresponds to the GLS behavior. However, when the effective range of interaction is 
decreased 7 > j c , the limiting case of local load sharing is approached and the strength 
of the system should vanish in the L — > 00 limit [2U [22]- In the isotropic case, 7 C falls 
in the vicinity of 7 C = 2 [16] . 

Figures [Ik., and [lb show the critical stress values ^ obtained for several ranges 
of interactions 7 and anisotropy strengths a. The data corresponds to systems with 
129 x 129 elements and a Weibull distribution of breaking thresholds with p = 2 and 



(a) (b) 




y 



Figure 1. The ultimate strength ^ of the system is studied, for a Weibull distribution 
of thresholds with p = 2 and ao = 1. a) We explore different ranges of interaction 7 
and anisotropy a. b) Results for different system sizes L and ranges of interaction 7 
are presented 
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o~ — 1. Two distinct regions can be easily identified in the graph [Th- For small 7, 
the critical strength — is independent, within statistical errors, of both the effective 

°~ o 

range of interaction 7 and the anisotropy strength a. As we have already pointed out, 
this behavior is expected for the standard GLS scheme. However, when the effective 
range of interaction decreases or the anisotropy strength increases, we found non-trivial 
dependencies. In Figure [Tb, the critical stress of a system with strong anisotropy 
(a ~ 0.9999) is detailed. We illustrate the model's behavior for several system sizes 
from L = 65 to L = 513. For small 7, — is independent of both effective range of 

O" o 

interaction and system size. However, as soon as the localized nature of the interaction 
becomes dominant, i.e. 7 > j c , ^ vanishes logarithmically with increasing system size. 
This qualifies for a genuine short-range behavior as found in local load sharing models, 
where the strength of the sample must vanish in the thermodynamic limit [241 [25] . In 
summary, in the regime 7 < 1, all the numerical findings are in excellent agreement 
with the mean field analytical prediction. Besides, figure [Th. suggests that the crossover 
point 7 C would differ from what was found using the isotropic stress transfer function 
[I~6] only for a very high anisotropy strength, i.e. a > 0.9999. 

Nevertheless, a priori one would expect the system size dependence to be more 
pronounced onces the anisotropy in introduced. Figure [lb shows that in the transition 
region the system size dependence is more appreciable than for the isotropic case [16]. On 
the left side of the transition region the critical stress ^ slowly increases with increasing 
system size. Note that for the isotropic case the convergence to the thermodynamic 
limit is faster; consequently, a more accurate estimation of the critical point 7 C = 2 
could be done [IB] . In the present case, in order to find a reasonable estimation of j c , 
we used the fact that on the short range interaction region (LLS) the convergence to 
the thermodynamic limit is qualitatively different than on the long range interaction 
region (GLS). On the GLS side of the transition region the critical stress ^ increases 
with increasing system size, towards the GLS exact solution. Contrary, on the LLS 
side of the transition region ^ vanishes logarithmically with increasing system size. 
Despite the fact that figure [EH seems to indicate 7 C 7^ 2 for a = 0.9999, figure [lb also 
suggests that for any value a < 1, the critical value 7 C shifts towards 7 C = 2 as the 
thermodynamic limit is reached. For elucidating if the value of 7 C changes onces the 
anisotropy is introduced we then focused on 7 = 1 and 7 = 2. Hence, changing the 
anisotropy strength a and studying the convergence of — to the thermodynamic limit, 

a 

a better estimation of 7 C is done. 

In Figure we describe the size dependency of the ultimate strength ^ for a 
Weibull distribution of thresholds with p = 2 and a = 1. Figure EH and figure 
[2b illustrate the outcomes at 7 = 1 and 7 = 2, respectively. Several anisotropy 
strengths were explored and it was found that when the system size is increased, the 
anisotropy plays a weaker role. Moreover, as expected the system size dependence 
is more pronounced for 7 = 2 than for 7 = 1. We note in Figj2b that numerical 
uncertainties surface when one gets numerically very close to a = 1, for small system 
sizes. That is due to the fact that in a 2D square lattice topology as a — > 1, the 
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load sharing from one row to the next might become so small that the rows would 
become effectively decoupled. That is certainly an undesirable topology effect which 
might also magnify the system size effects. However, our results indicate that even in 
the presence of high anisotropy, at 7 C = 2, the system shows a tendency to behave as 
GLS when approaching to the thermodynamic limit, within our numerical uncertainties. 
Consequently, the crossover point for the anisotropic variable range of interaction given 
by Eq.(j3]) results at 7 C = 2. 

To test the universality of this statement, we did calculations for several threshold 
distributions. A Weibull distribution with p = 10, as well as a uniform distribution were 
used for comparison. Moreover, to access accurately the crossover point, several system 
sizes were considered. In every case, we changed a using the values (0.9; 0.99; 0.999 and 
0.9999). Our aim is to elucidate how far from GLS behavior the system is for 7 = 2, 
after the anisotropy is introduced. As we pointed out before, this value 7 C = 2 defines 
the crossover between GLS and LLS behaviors, for the isotropic case. 

In Figure El results of the ultimate strength of samples with strong anisotropy are 
shown in detail. The scaled magnitude x = (a c — a)L^ is plotted against the "distance" 
5a = ctgls — fr° m the well-defined GLS's behavior. Note that each symbol (i.e. 
square, circle and diamond) is related to a different threshold distribution. In addition, 
the sizes of the symbols are linked to different system sizes (the bigger the symbol, the 
larger the system size). Different regions of the plot, illustrated with the same symbol in 
five different sizes, correspond to a given value of a (from left to the right 0.9; 0.99; 0.999 
and 0.9999). In the plots the data, corresponding to each threshold distribution, are 
aligned finding data collapse, 

5a = (a GLS - ~ F((a c - a)L*) ~ F(x) (6) 

where we have introduced the scaling function F((a c — a)L^), with £ = 0.5 and 




Figure 2. We show results of ultimate strength as a function of the system size, for 
several anisotropy strengths. A Weibull distribution of thresholds with p = 2 and 
(To = 1 is used, a) results obtained at 7 = 1, b) results obtained at 7 = 2 
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Figure 3. The ultimate strength of the system o c is studied at 7 = 2, through 
the magnitude 8a — cfgls — —■ Several anisotropy strength values and [L = 
65, 129, 257, 513, 1025) are illustrated. 



a c — 1. The results for a uniform distribution and Weibull with p = 2 are very similar. 
Furthermore, decreasing the disorder, e.g. Weibull with p = 10 and a Q = 1, only 
magnifies the topology and finite size effects. The scaling exponent £ ~ 0.5 is universal, 
defining a new crossover exponent in a and L. From those results, one can conclude 
that the anisotropy does not change the crossover point for the range of interaction in 
2D, 7 C = 2. For finite systems, the global behavior is reached very slowly, as we get far 
from a c = 1. We also notice that the distance to the GLS behavior 5a = {o~gls ~ f 2 ) 
does vanish in the infinite system size limit, as a power law 5a ~ , and we estimate 
f3 = 0.37 ± 0.03. We conclude, that even in presence of high anisotropy the behavior of 
the system for 7 < 7 C = 2 shows a tendency to GLS as L — > 00. The infinite system 
would display a qualitatively different behavior only for a c = 1, where it would become 
effectively a ID model, with j c — 1. 

The fracture process can also be described by the precursory activity before 
complete breakdown. The statistical properties of rupture sequences are characterized 
by the avalanche size distribution. From an experimental point of view the precursory 
activity is related to the acoustic emissions generated during the fracture of materials 
[2T| I2"SI |2"§1 [3D] . The avalanche size distribution is a measure of causally connected 
broken sites. All the intact elements have a non-zero chance to fail independently of 
the (spatial) rupture history, and any given element could be near to its rupture point 
regardless of its position in the lattice. 
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Figure 4. Avalanche size distribution for N = 257 x 257 fibers. Results corresponding 
to several anisotropic strengths are illustrated 



The highly fluctuating activity is certainly related to the long range interactions 
where the avalanche size distribution can usually be well fitted by a power law 
.P(A) ~ A - 2. This actually corresponds to the mean field scenario, GLS [22l 1231 l24"t [25] . 
However, when the spatial correlations are important LLS, stress concentration takes 
place in the elements located at the perimeter of an already formed cluster. Hence, 
elements far away from the clusters of broken elements have significantly lower stresses 
and thus the size of the largest avalanche is reduced as well as the number of failed 
elements belonging to the same avalanche, leading to lower precursory activity. 

Figure H] illustrates the avalanche statistics obtained for systems with N = 257x257 
fibers. In every case, we have set 7 = 2 and used different values of the anisotropy 
strength a. It is noticeable that the avalanche size distributions can always be fitted 
to a power law with a non-trivial exponent. However, as we get far from a = 1 the 
exponent tends asymptotically to the value r = §, reflecting the tendency to recover 
the GLS's behavior. 

To characterize the system size dependence, we propose the following scaling ansatz 
for the avalanche size distribution, 

D(A,L) = A C 7 0(A) (7) 

where A^ is a characteristic avalanche A c h ~ L e and g(x) is a scaling function that goes 
like g (x) = x~ T for x < 1 and decays faster than a power law for x > 1 . In Figure \5[ the 
avalanche statistics resulting from systems with different sizes are shown. The model is 
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Figure 5. The normalized distributions D(A) of avalanches for a = 0.999 are 
illustrated, several system sizes are considered. In the inset we verify the scaling 
ansatz given by Eq.([7]) with A c /j ~ L e . 

again investigated at 7 = 2, and very strong anisotropy a = 0.999. The scaled function 
is presented in the inset. It can be seen that the data collapse yields 6 = 0.6 ±0.1 and 
r = 2.50 ± 0.02 with a power law over several orders of magnitude in This result 
suggests that for 7 = 2, even in presence of very high anisotropy, the avalanches are 
distributed as in the case of GLS, namely as D(A) ~ A~^. Thus, the crossover point is 
still at 7 C = 2, even in the case of an anisotropic stress redistribution (given by Eq.Q). 

4. Discussion 

Long-range fiber bundle models can be considered as a first approximation to model 
the fracture behavior of a disordered elastic medium. It has been proposed in several 
instances to replace the full solution of the elastic equations by a Green function [311 [32] . 
This method has the advantage to avoid the computational cost involved in the inversion 
of the elastic equations, but in principle it is only accurate for diluted damage. A 
particularly simple example is provided by the random fuse model (RFM) in which 
a lattice of conducting bonds with random failure thresholds is loaded by applying 
an external voltage at the two ends of the lattice. When a fuse fails the current is 
redistributed to the neighboring fuses by solving the Kirchhoff equations. When only a 
few bonds are broken, the current is transferred according to the homogeneous lattice 
Green function which is given by F(r) = x/r 3 in two dimensions [31]. It is therefore 
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Figure 6. Global strength distribution of samples with several types of load sharing 
(L = 129, averaged over 5000 different configurations). GLS and RFM are very the 
similar. 



a long-range (with 7 = 2) and anisotropic load transfer function, but of a different 
character than the one we have used here. 

We have simulated a fiber bundle model using a load transfer function inspired 
by the RFM. The general result is that the failure properties of our model resemble 
very much those of GLS fiber bundles rather than those of the original RFM. This is 
particularly apparent when looking at the strength distribution, displayed in Figure El 
Both the mean-field GLS approach and our result obtained with the RFM load transfer 
function obey Gaussian statistics. Notice that the original RFM displays instead a 
qualitatively different log- normal strength distribution [3J. This confirms that the Green 
function approach is reliable at most in the initial stages of the damage accumulation 
process and it is not correct to describe the global failure of the RFM. Figure [6] also 
shows that anisotropic LLS functions result instead in a larger scattering of ultimate 
strength, and their asymmetric distributions are usually fitted by Weibull distribution 
functions. It is noticeable from the data that introducing very high anisotropy strength 
amplifies the scattering in the global strength of the samples. However, the strength 
distribution for the anisotropic system appears to be closer to a Gaussian in contrast 
to the classical Weibull behavior, which is usually obtained for LLS approaches. 

In conclusion, we have studied a discrete fracture model where the interaction 
among elements is considered to decay anisotropically with the distance from an intact 
element to the rupture point. The two classical regimes (local and global) are found 
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as the exponent of the stress-transfer function varies and a crossover point is again 
identified in the vicinity of 7 C = 2. The strength of the material for 7 < 7 C does not 
depend on both the system size and 7 qualifying for mean-field behavior, whereas for 
the short range regime, the critical load vanishes in the thermodynamic limit. The 
behavior of the model at both sides of the crossover point was numerically studied 
by recording the avalanche and the critical stress for several system sizes. From our 
numerical results, one can certainly conclude that the anisotropy does not change the 
crossover point 7 C = 2 in the 2D model, in the infinite system size limit. The 2D model 
would display a qualitatively different behavior only for a c — 1, where it would become 
effectively a ID model, with 7 C = 1. Moreover, in finite systems for 7 < 2, the global 
load sharing behavior is very slowly recovered as we get far from a c — 1, within our 
numerical uncertainties. 
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